//*******************************************************************
// For w control
//*******************************************************************
#include "plant_1DTRI_RDM.h"

const double TOLERANCE = 1e-15;
const double VERSION   = 1.1;

//*******************************************************************
// Constant defined for debug
//*******************************************************************
const double T_START	= 0.0;

//*******************************************************************
//	PUBLIC FUNCTIONS
//  - Constructor/Destructor
//  - Simulation control
//    * reset
//    * update
//    * notfinish
//*******************************************************************
// must reset() before running the kMC model
lattice_1DTRI::lattice_1DTRI(int iLatticeSize, double iTimeEnd, double iTimeInt, double iTempIni, double iDepRateIni, int iMIni, int iMInc, char *iFileName, int iDetail)
{
	int succ=1;
	int i;
	temp_ini = iTempIni;
	dep_rate_ini = iDepRateIni;
	t_start = T_START;
	t_final = iTimeEnd;
	t_int = iTimeInt;
	LatticeSize = iLatticeSize;
	MLyr=iMIni;
	M_ini = iMIni;
	M_inc = iMInc;
		
	hts			= new int[LatticeSize];
	hts_half	= new int[LatticeSize];
	occu		= new int[LatticeSize];
	phi_nx		= new double[MODE*2*LatticeSize];
	psi_nx		= new double[MODE*2*LatticeSize];

	record.FileName = iFileName;
	record.rec_num	= int((t_final-t_start)/t_int)+1;
	record.rec_idx	= record.rec_num;	// this assignment guarantees notfinished() returns false
	record.pRecords = new struct OneRecord[record.rec_num];

	init_phi_psi();
	init_prev_next();
	initiate_stats();

	if(!succ)
	{
		cerr << "lattice_1DTRI::Initialization fails\n";
	}
	
	Detail = iDetail;
	if(Detail==1)
	{
		fd_time.open("KMC_time_d.dat");
		fd_T.open("KMC_T_d.dat");
		fd_W.open("KMC_W_d.dat");
		fd_SOR.open("KMC_SOR_d.dat");	
		fd_H.open("KMC_H_d.dat");
		fd_R2.open("KMC_R2_d.dat");
		fd_M2.open("KMC_M2_d.dat");
		//fd_alpha.open("KMC_alpha_d.dat");
		//fd_beta.open("KMC_beta_d.dat");
	}
	srand((unsigned)time( NULL ) );
	return;
}
//*******************************************************************
lattice_1DTRI::~lattice_1DTRI()
{
	avg_stats();
	// -------------------------------
	ofstream fout;
	int i;

	fout.open(record.FileName);
	fout << "class: lattice_1DTRI\n";
	fout << "LatticeSize = " << LatticeSize << "\n";
	fout << "TopLyr      = " << TopLyr << "\n";
	fout << "Temperature = " << Temperature << "\n";
	fout << "DIFF_E0     = " << DIFF_E0 << "\n";
	fout << "DIFF_R0     = " << DIFF_R0 << "\n";
	fout << "Repeated    = " << Repeated << "\n";
	fout << "T_START     = " << t_start << "\n";
	fout << "T_END       = " << t_final << "\n";
	fout << "T_INT       = " << t_int << "\n";
	fout << "T_NODE      = " << record.rec_num << "\n";
	fout << "M_ini       = " << M_ini << "\n";
	fout << "M_inc       = " << M_inc << "\n";
	//fout << setw(6) << "t(s)" << setw(15) << "<T>(K)" << setw(15) << "<W>(layer/s)" << setw(15) << "<rho_s>" << setw(15) << "<rho_s^2>";
	fout << setw(6) << "t(s)" << setw(15) << "<T>(K)" << setw(15) << "<W>(layer/s)";
	fout << setw(15) << "<rho_g>" << setw(15) << "var(rho_g)" << setw(15) << "<h>" << setw(15) << "var(h)" << setw(15) << "<R2>" << setw(15) << "var(R2)" << setw(15) << "<slope>" << setw(15) << "var(slope)" << "\n";
	for(i = 0;i<record.rec_num;i++)
	{
		fout << setw(6)  << setprecision(6)<< record.pRecords[i].time;
		fout << setw(15) << setprecision(6)<< record.pRecords[i].T;
		fout << setw(15) << setprecision(6)<< record.pRecords[i].w;
		fout << setw(15) << setprecision(6)<< record.pRecords[i].SOR_mean;
		fout << setw(15) << setprecision(6)<< record.pRecords[i].SOR_var;
		fout << setw(15) << setprecision(6)<< record.pRecords[i].H_mean;
		fout << setw(15) << setprecision(6)<< record.pRecords[i].H_var;
		fout << setw(15) << setprecision(6)<< record.pRecords[i].R2_mean;
		fout << setw(15) << setprecision(6)<< record.pRecords[i].R2_var;
		fout << setw(15) << setprecision(6)<< record.pRecords[i].M2_mean;
		fout << setw(15) << setprecision(6)<< record.pRecords[i].M2_var;
		if(i<record.rec_num-1) {fout << "\n";}
	}
	fout.close();
	if(Detail==1) 
	{
		fd_time.close();
		fd_T.close();
		fd_W.close();
		fd_SOR.close();
		fd_H.close();
		fd_R2.close();
		fd_M2.close();
		//fd_alpha.close();
		//fd_beta.close();
	}
	// -------------------------------
	fin_lattice();
	
	// Dimension = LatticeSize
	delete [] PREV_ST;
	delete [] NEXT_ST;
	delete [] hts;
	delete [] hts_half;
	delete [] occu;
	// Dimension = MODE*(2*LatticeSize)
	delete [] phi_nx;
	delete [] psi_nx;
	// Dimension = record.rec_num
	delete [] record.pRecords;

	return;
}
//*******************************************************************
void lattice_1DTRI::reset()
{
	if(notfinish()) {cerr << "ERROR@lattice_1DTRI::reset() No resetting when the current kMC simulation not finished. rec_idx = " << record.rec_idx << "\n";system("pause");exit(0);}

	fin_lattice();

	Repeated++;
	
	init_lattice();
	Dep_rate = dep_rate_ini;
	Temperature = temp_ini;	// always init_diff_rate() and init_p_diff() after changing Temperature
	init_diff_rate();
	init_p_diff();
	record.rec_idx = 0;
	update_records(0.0);	// update dep_time and all records at time 0
	return;
}
//*******************************************************************
void lattice_1DTRI::get()
{
	cout << "class: lattice_1DTRI\n";
	cout << "MLyr        = " << MLyr << "\n";
	cout << "TopLyr      = " << TopLyr << "\n";
	cout << "LatticeSize = " << LatticeSize << "\n";
	cout << "M_ini       = " << M_ini << "\n";
	cout << "M_inc       = " << M_inc << "\n";
	cout << "record.rec_num      = " << record.rec_num << "\n";
	cout << "dep_time    = " << dep_time << "\n";

	return;
}
//*******************************************************************
double lattice_1DTRI::get_dep_time()
{
	return dep_time;
}
//*******************************************************************
struct lattice_state lattice_1DTRI::get_lattice_state()
{
	return LState;
}
//*******************************************************************
// 2009.07.29 Gangshi: a new algorithm which will move a fixed time step, whether a kMC event is included or not
// 2009.07.30 Gangshi: introduce DepDuration as the duration of DepRate as a constant input to the kMC simulation
void lattice_1DTRI::update(double iTemperature, double iDepRate, double iDepDuration)
{
	double rdm;
	double rate_ttl;
	double dep_time_next;
	double t_MC;

	if(iDepRate < 0.0 || iDepRate > 1000) {cerr << "ERROR@lattice_1DTRI::update() Illegal value of iDepRate = " << iDepRate << "\n";system("pause");exit(0);}
	if(iTemperature < 300.0 || iTemperature > 900.0) {cerr << "ERROR@lattice_1DTRI::update() Illegal value of iTemperature = " << iTemperature << "\n";system("pause");exit(0);}

	Dep_rate = iDepRate;
	if(Temperature!=iTemperature)
	{
		Temperature = iTemperature;	// always init_diff_rate() and init_p_diff() after changing Temperature
		init_diff_rate();
		init_p_diff();
	}
	
	// repeat kMC steps and update state vectors when possible until dep_time>dep_time_next, the next should-be dep_time
	// when going beyond, quit while loop without running the kMC event, since no event happened in this loop (the time of next event is beyond time_next)
	dep_time_next=dep_time+iDepDuration;
	while(1)
	{
		rate_ttl = Dep_rate*LatticeSize+diff_ttl;
		t_MC=OneMCTime(rate_ttl);
		if(rate_ttl<0.0) {cerr << "ERROR@lattice_1DTRI::update() Illegal (negative) value of rate_ttl = " << rate_ttl << "\n";system("pause");exit(0);}
		if (rate_ttl==0.0 || dep_time+t_MC>dep_time_next) break;
		// run one kMC step
		rdm=RandRand()*rate_ttl;
		if(rdm<=diff_ttl)
		{
			diffusion();
		}
		else
		{
			deposition();
		}
		update_records(dep_time+t_MC);	// update records each time dep_time is changed
	}
	try
	{
		update_records(dep_time_next);	// update records each time dep_time is changed
	}
	catch(const char* s)
	{
		cerr << "\n" << s;
		cerr << "\nError when calling: update_records(dep_time_next)";
		system("pause");
	}

	LState.w = iDepRate;
	LState.top = TopLyr;
	return;
}
//*******************************************************************
void lattice_1DTRI::update2(double time, double iW,double iT)
{
	if(time<dep_time)
		return;
	else
	{
		update(iT,iW,time-dep_time);
	}
}

//*******************************************************************
void lattice_1DTRI::display_lattice()
{
	int i,j;
	printf("----------------------------------------------------------\n");
	for(j=MLyr+2;j>=0;j--)
	{
		for(i=0;i<LatticeSize;i++)
		{
			if(j&0x1)
			{
				printf(" %s",pp_lattice[j][i]&0x1?"X":"0");
			}
			else
			{
				printf("%s ",pp_lattice[j][i]&0x1?"X":"0");
			}
		}
		printf("\n");
	}
	return;
}
//*******************************************************************
int lattice_1DTRI::notfinish()
{
	return (record.rec_idx < record.rec_num);	
}
//*******************************************************************

//*******************************************************************
//
//	PRIVATE FUNCTIONS
//
//  - Initialization
//  - Lattice operation (Level 1)
//    * add_particle
//    * remove_particle
//    * add_relax()
//    * get_relative_sites
//    * get_relative_site
//  - Lattice operation (Level 2)
//    * deposition
//    * diffusion
//  - Utilities
//    * Random number generator
//  - I/O routine
//*******************************************************************

//*******************************************************************
// Private functions/Initialization
//*******************************************************************
void lattice_1DTRI::init_lattice()
{
	int succ;
	int i,j;

	MLyr = M_ini;

	// Initialize Lattice, modified from void init_lattice()
	pp_lattice=(unsigned int **)malloc((M_ini+3)*sizeof(unsigned int *));
	if(pp_lattice==NULL)
		succ=0;
	else
	{
		pp_lattice[0]=(unsigned int *)malloc((M_ini+3)*LatticeSize*sizeof(unsigned int));
		if(pp_lattice[0]==NULL)
			succ=0;
		else
		{
			memset(pp_lattice[0],0,(M_ini+3)*LatticeSize*sizeof(unsigned int));
			for(j=1;j<M_ini+3;j++)
			{
				pp_lattice[j]=pp_lattice[0]+j*LatticeSize;
			}
			for(i=0;i<LatticeSize;i++)
			{
				// initialize the fixed layer
				pp_lattice[0][i]=ST_FIX;
				// initialize the sites above the fixed layer
				pp_lattice[1][i]=0x0000000c;
			}
			TopLyr=0;
			//initialize the diffusion rate table
			p_diff=(double *)malloc((M_ini+3)*sizeof(double));
			if(p_diff==NULL)
				succ=0;
			else
				memset(p_diff,0,(M_ini+3)*sizeof(double));
		}
	}

	return;
}
//*******************************************************************
void lattice_1DTRI::init_prev_next()
{
	int i;

	PREV_ST = new int[LatticeSize];
	PREV_ST[0]=LatticeSize-1;
	for(i=1;i<LatticeSize;i++)
	{
		PREV_ST[i]=i-1;
	}

	NEXT_ST = new int[LatticeSize];
	NEXT_ST[LatticeSize-1]=0;
	for(i=0;i<LatticeSize-1;i++)
	{
		NEXT_ST[i]=i+1;
	}

	return;
}
//*******************************************************************
void lattice_1DTRI::init_p_diff()
{
	int i,j;

	diff_ttl=0.0;
	for(j=1;j<MLyr+1;j++)
	{
		p_diff[j]=0.0;
		for(i=0;i<LatticeSize;i++)
		{
			if(GET_SITE(j,i,0))
			{
				p_diff[j]+=DIFF_RT[GET_NB_NUM(j,i)];
			}
		}
		diff_ttl+=p_diff[j];
	}

	return;
}
//*******************************************************************
void lattice_1DTRI::init_diff_rate()
{
	int i;

	DIFF_RT[0]=0;
	for(i=1;i<6;i++)
	{
		DIFF_RT[i]=DIFF_R0*exp(-i*(DIFF_E0+0.0)*EV2KB/(Temperature+0.0));
	}
	DIFF_RT[6]=0;

	return;
}
//*******************************************************************
void lattice_1DTRI::init_phi_psi()
{
	int i,j;
	double n,x;
	
	for (i=0;i<MODE;i++)
	{
		n=i+1.0;
		for (j=0;j<2*LatticeSize;j++)
		{
			x=-PI+j*PI/LatticeSize;
			phi_nx[i*2*LatticeSize+j] = RECP_SQRT_PI*sin(n*x);
			psi_nx[i*2*LatticeSize+j] = RECP_SQRT_PI*cos(n*x);
		}
	}

	return;
}
//*******************************************************************
//initialize states and variables of the lattice before any loop runs
void lattice_1DTRI::initiate_stats()
{
	int i;

	cout << "\nt_int = " << t_int;
	for (i=0;i<record.rec_num;i++)
	{
		record.pRecords[i].time			= t_start+i*t_int;
		record.pRecords[i].T			= 0.0;
		record.pRecords[i].w			= 0.0;
		record.pRecords[i].H_mean		= 0.0;
		record.pRecords[i].H_var		= 0.0;
		record.pRecords[i].R2_mean		= 0.0;
		record.pRecords[i].R2_var		= 0.0;
		record.pRecords[i].SOR_mean		= 0.0;
		record.pRecords[i].SOR_var		= 0.0;
		record.pRecords[i].pSOR_mean	= 0.0;
		record.pRecords[i].pSOR_var		= 0.0;
		record.pRecords[i].M2_mean	= 0.0;
		record.pRecords[i].M2_var	= 0.0;
	}
	// initialize other variables
	Repeated=0;
	dep_time=-1;
	pp_lattice=NULL;
	p_diff=NULL;

	return;
}

//*******************************************************************
// Private functions/Lattice operation (level 1)
//*******************************************************************
void lattice_1DTRI::fin_lattice()
{
	int j;
	// Dimension with LatticeSize*MLyr
	if(pp_lattice!=NULL)
	{
		for(j=M_INI+3;j<MLyr+3;j+=M_INC)
		{
			if (pp_lattice[j]!=NULL)
			{
				free(pp_lattice[j]);
				pp_lattice[j]=NULL;
			}
		}
		if (pp_lattice[0]!=NULL)
		{
			free(pp_lattice[0]);
			pp_lattice[0]=NULL;
		}
		free(pp_lattice);
		pp_lattice=NULL;
	}
	
	if(p_diff!=NULL)
	{
		free(p_diff);
		p_diff=NULL;
	}

	return;
}
//*******************************************************************
int lattice_1DTRI::inc_lattice()
{
	int j,MLyr_old;
	int succ=1;
	unsigned int	**pp_lattice_tmp;
	double	*p_diff_tmp;

	MLyr_old=MLyr;
	MLyr+=M_INC;

	pp_lattice_tmp=(unsigned int **)malloc((MLyr+3)*sizeof(unsigned int *));
	if (pp_lattice_tmp == NULL)
		succ=0;
	else
	{
		pp_lattice_tmp[MLyr_old+3]=(unsigned int *)malloc(M_INC*LatticeSize*sizeof(unsigned int));
		if(pp_lattice_tmp[MLyr_old+3]==NULL)
		{
			free(pp_lattice_tmp);
			succ=0;
		}
		else
		{
			memset(pp_lattice_tmp[MLyr_old+3],0,M_INC*LatticeSize*sizeof(unsigned int));
			for(j=0;j<M_INC;j++)
			{
				pp_lattice_tmp[MLyr_old+3+j]=pp_lattice_tmp[MLyr_old+3]+j*LatticeSize;
			}
			memcpy(pp_lattice_tmp,pp_lattice,(MLyr_old+3)*sizeof(unsigned int *));
			free(pp_lattice);
			pp_lattice=pp_lattice_tmp;
			pp_lattice_tmp=NULL;

			//increase the diffusion rate table
			p_diff_tmp=(double *)malloc((MLyr+3)*sizeof(double));
			if(p_diff==NULL)
				succ=0;
			else
				memset(p_diff_tmp,0,(MLyr+3)*sizeof(double));
			memcpy(p_diff_tmp,p_diff,(MLyr_old+3)*sizeof(double));
			free(p_diff);
			p_diff=p_diff_tmp;
			p_diff_tmp=NULL;
		}
	}

	return succ;
}
//*******************************************************************
void lattice_1DTRI::get_relative_sites(int lyr, int st, int *nb_lyrs, int *nb_sts)
{
	nb_lyrs[0]=lyr;			nb_sts[0]=st;
	if (lyr&0x1)
	{
		nb_lyrs[1]=lyr;		nb_sts[1]=PREV_ST[st];
		nb_lyrs[2]=lyr-1;	nb_sts[2]=st;
		nb_lyrs[3]=lyr-1;	nb_sts[3]=NEXT_ST[st];
		nb_lyrs[4]=lyr;		nb_sts[4]=NEXT_ST[st];
		nb_lyrs[5]=lyr+1;	nb_sts[5]=NEXT_ST[st];
		nb_lyrs[6]=lyr+1;	nb_sts[6]=st;
	}
	else
	{
		nb_lyrs[1]=lyr;		nb_sts[1]=PREV_ST[st];
		nb_lyrs[2]=lyr-1;	nb_sts[2]=PREV_ST[st];
		nb_lyrs[3]=lyr-1;	nb_sts[3]=st;
		nb_lyrs[4]=lyr;		nb_sts[4]=NEXT_ST[st];
		nb_lyrs[5]=lyr+1;	nb_sts[5]=st;
		nb_lyrs[6]=lyr+1;	nb_sts[6]=PREV_ST[st];
	}

	return;
}
//*******************************************************************
void lattice_1DTRI::get_relative_site(int lyr, int st, int nb, int *nb_lyr, int *nb_st)
{
	switch(nb)
	{
		case 1: case 4: *nb_lyr=lyr; break;
		case 2: case 3: *nb_lyr=lyr-1; break;
		case 5: case 6: *nb_lyr=lyr+1; break;
	}
	if (lyr&0x1)
	{
		switch(nb)
		{
			case 1: *nb_st=PREV_ST[st]; break;
			case 2: case 6: *nb_st=st; break;
			case 3: case 4: case 5: *nb_st=NEXT_ST[st]; break;
		}
	}
	else
	{
		switch(nb)
		{
			case 1: case 2: case 6: *nb_st=PREV_ST[st]; break;
			case 3: case 5: *nb_st=st; break;
			case 4: *nb_st=NEXT_ST[st]; break;
		}
	}
	return;
}
//*******************************************************************
void lattice_1DTRI::add_particle(int lyr, int st)
{
	int n;
	int lyrs[7], sts[7];
	int nb_num;

	TopLyr=TopLyr>lyr?TopLyr:lyr;
	while(TopLyr>MLyr) inc_lattice();

	get_relative_sites(lyr, st, lyrs, sts);

	// test if the site is vacant and give warning if occupied
	if(GET_SITE(lyrs[0],sts[0],0))
	{
		printf("the site (%d, %d) is occupied and no particles can be added to the site!\n",lyrs[0],sts[0]);
	}

	SET_SITE(lyrs[0],sts[0],0);
	for(n=1;n<4;n++)
	{
		SET_SITE(lyrs[n],sts[n],n+3);
	}
	for(n=4;n<7;n++)
	{
		SET_SITE(lyrs[n],sts[n],n-3);
	}

	// remove the sum of diffusion rate at lyr-1, lyr and lyr+1 layer
	// if lyr is the bottom layer, no removal
	diff_ttl-=p_diff[lyr+1];
	diff_ttl-=p_diff[lyr];
	if(lyr>1)
	{
		diff_ttl-=p_diff[lyr-1];
	}
	// Adjust diffusion rates of the relative sites
	// if the relative site is on the fixed layer, no adjustment
	for(n=1;n<7;n++)
	{
		if(GET_SITE(lyrs[n],sts[n],0)&&(lyrs[n]!=0))
		{
			nb_num=GET_NB_NUM(lyrs[n],sts[n]);
			p_diff[lyrs[n]]-=DIFF_RT[nb_num-1];
			p_diff[lyrs[n]]+=DIFF_RT[nb_num];
		}
	}
	p_diff[lyr]+=DIFF_RT[GET_NB_NUM(lyr,st)];
	// add the sum of diffusion rate at lyr-1, lyr and lyr+1 layer
	// if lyr is the bottom layer, no addtion
	diff_ttl+=p_diff[lyr+1];
	diff_ttl+=p_diff[lyr];
	if(lyr>1)
	{
		diff_ttl+=p_diff[lyr-1];
	}

	return;
}
//*****************************************************************************
void lattice_1DTRI::remove_particle(int lyr, int st)
{
	int n;
	int lyrs[7], sts[7];
	int nb_num;

	get_relative_sites(lyr, st, lyrs, sts);

	// test if the site is occupied and give warning if vacant
	if(!GET_SITE(lyrs[0],sts[0],0))
	{
		printf("the site (%d, %d) is vacant and no particles can be removed to the site!\n",lyrs[0],sts[0]);
	}

	RESET_SITE(lyrs[0],sts[0],0);
	for(n=1;n<4;n++)
	{
		RESET_SITE(lyrs[n],sts[n],n+3);
	}
	for(n=4;n<7;n++)
	{
		RESET_SITE(lyrs[n],sts[n],n-3);
	}

	// remove the sum of diffusion rate at lyr-1, lyr and lyr+1 layer
	// if lyr is the bottom layer, no removal
	diff_ttl-=p_diff[lyr+1];
	diff_ttl-=p_diff[lyr];
	if(lyr>1)
		diff_ttl-=p_diff[lyr-1];
	// Adjust diffusion rates of the relative sites
	// if the relative site is on the fixed layer, no adjustment
	for(n=1;n<7;n++)
	{
		if(GET_SITE(lyrs[n],sts[n],0)&&(lyrs[n]!=0))
		{
			nb_num=GET_NB_NUM(lyrs[n],sts[n]);
			p_diff[lyrs[n]]-=DIFF_RT[nb_num+1];
			p_diff[lyrs[n]]+=DIFF_RT[nb_num];
		}
	}
	p_diff[lyr]-=DIFF_RT[GET_NB_NUM(lyr,st)];
	// add the sum of diffusion rate at lyr-1, lyr and lyr+1 layer
	// if lyr is the bottom layer, not addition
	diff_ttl+=p_diff[lyr+1];
	diff_ttl+=p_diff[lyr];
	if(lyr>1)
		diff_ttl+=p_diff[lyr-1];

	return;
}
//*******************************************************************
int lattice_1DTRI::add_relax(int lyr, int st)
{
	int i;
	int nb_num;
	int lyrs[7],sts[7],nb_nums[7];
	int max_nb,max_nb_num,max_lyrs[5],max_sts[5];
	int lyr_relax,st_relax;
	int rdm;

	nb_num=GET_NB_NUM(lyr,st);

	if(GET_SITE(lyr,st,0))
		return 0;

	if(nb_num==0)
		return 0;

	if(!RELAX)
	{
		add_particle(lyr,st);
		return 1;
	}
	// If a site has only one neighbor, then find the neighboring vacant with maximum neighboors
	// If there are multiple neighboring vacant with the same amount of neighbors, then select 
	// one randomly
	if(nb_num==1)
	{
		max_nb=2;
		get_relative_sites(lyr, st, lyrs, sts);

		for(i=1;i<7;i++)
		{
			if(!GET_SITE(lyrs[i],sts[i],0))
			{
				nb_nums[i]=GET_NB_NUM(lyrs[i],sts[i]);
				if(nb_nums[i]>max_nb)
					max_nb=nb_nums[i];
			}
			else
			{
				nb_nums[i]=-1;
			}
		}
		max_nb_num=0;
		for(i=1;i<7;i++)
		{
			if(nb_nums[i]==max_nb)
			{
				max_lyrs[max_nb_num]=lyrs[i];
				max_sts[max_nb_num]=sts[i];
				max_nb_num++;
			}
		}
		rdm=(int)(RandRand()*max_nb_num);
		lyr_relax=max_lyrs[rdm];
		st_relax=max_sts[rdm];
	}
	else
	{
		lyr_relax=lyr;
		st_relax=st;
	}
	

	add_particle(lyr_relax,st_relax);


	return 1;
}

//*******************************************************************



//*******************************************************************
// Private functions/Lattice Operation (level 2)
//*******************************************************************
void lattice_1DTRI::deposition()
{
	double x0,ang;
	double sinang,cosang;
	double coll_width;
	double x_coll,x_colls[3];
	int coll_lyr_num=0,test_lyr_num=0;
	int lyr,st,st_left,st_right,st_stat;
	int coll_lyrs[3],coll_sts[3];
	int coll_lyr,coll_st,coll_nb;
	int dep_lyr,dep_st;
	double coll_dist,coll_dists[3];
	double dx,dy;

	int st_lefts[3],st_rights[3];
	unsigned int st_stats[3];

	x0=RandRand()*LatticeSize;
	//ang=(RandRand()-0.5)*PI;
	ang=0.0/180.0*PI; // incident angle 0 deg // verticle incidence

	sinang=sin(ang);
	cosang=cos(ang);

	coll_width=1.0/cosang;

	lyr=TopLyr;
	while((test_lyr_num<3)&&(lyr!=-1))
	{
		x_coll=x0+(MLyr+2-lyr)*H_ML*sinang/cosang;
		st_left=(int)(ceil(x_coll-0.5*(lyr&0x1)-coll_width));
		st_right=(int)(floor(x_coll-0.5*(lyr&0x1)+coll_width));
		if(ang>0)
		{
			st=st_left;
			while(!(st_stat=GET_SITE(lyr,(st%LatticeSize+LatticeSize)%LatticeSize,0))&&(st!=st_right)) st=st++;
		}
		else
		{
			st=st_right;
			while(!(st_stat=GET_SITE(lyr,(st%LatticeSize+LatticeSize)%LatticeSize,0))&&(st!=st_left)) st=st--;
		}
		if(st_stat&&(GET_NB_NUM(lyr,(st%LatticeSize+LatticeSize)%LatticeSize)!=6))
		{
			coll_lyrs[coll_lyr_num]=lyr;
			coll_sts[coll_lyr_num]=st;
			x_colls[coll_lyr_num]=x_coll;
			st_lefts[coll_lyr_num]=st_left;
			st_rights[coll_lyr_num]=st_right;
			st_stats[coll_lyr_num]=pp_lattice[lyr][(st%LatticeSize+LatticeSize)%LatticeSize];
			coll_lyr_num++;
			test_lyr_num++;
		}
		else if(coll_lyr_num)
			test_lyr_num++;
		lyr--;
	}

	if(coll_lyr_num<=0)
	{
		return;
	}

	for(lyr=0;lyr<coll_lyr_num;lyr++)
	{
		dx=x_colls[lyr]-(coll_sts[lyr]+0.5*(coll_lyrs[lyr]&0x1));
		coll_dists[lyr]=(MLyr+2-coll_lyrs[lyr])*H_ML/cosang-dx*sinang-sqrt(1-dx*dx*cosang*cosang);
	}
	coll_lyr=coll_lyrs[0];
	coll_st=coll_sts[0];
	coll_dist=coll_dists[0];
	for(lyr=1;lyr<coll_lyr_num;lyr++)
		if(coll_dists[lyr]<coll_dist)
		{
			coll_lyr=coll_lyrs[lyr];
			coll_st=coll_sts[lyr];
			coll_dist=coll_dists[lyr];
		}
	dx=(x0+coll_dist*sinang)-(coll_st+0.5*(coll_lyr&0x1));
	dy=((MLyr+2)*H_ML-coll_dist*cosang)-coll_lyr*H_ML;
	if(dy>0.5)
	{
		if(dx>0)
			coll_nb=5;
		else if(dx<0)
			coll_nb=6;
		else if(ang>0)
			coll_nb=5;
		else
			coll_nb=6;
	}
	else if((dy<=0.5)&&(dy>=-0.5))
	{
		if(dx>0)
			coll_nb=4;
		else
			coll_nb=1;
	}
	else
	{
		if(dx>0)
			coll_nb=3;
		else
			coll_nb=2;
	}

	coll_st=coll_st%LatticeSize;
	if(coll_st<0) coll_st+=LatticeSize;
	if(GET_SITE(coll_lyr,coll_st,coll_nb))
	{
		if(coll_nb==4)
			coll_nb=5;
		else if(coll_nb==1)
			coll_nb=6;
		else
			printf("Impossible collision!\n");
	}

	get_relative_site(coll_lyr,coll_st,coll_nb,&dep_lyr,&dep_st);

	if(GET_SITE(dep_lyr,dep_st,0))
	{
		printf("Cannot deposit onto the neighbor site of (%d,%d)!\n",dep_lyr,dep_st);
	}

	add_relax(dep_lyr,dep_st);
	return;
}

//*******************************************************************
int lattice_1DTRI::diffusion()
{
	int lyr,st,n;
	int nb_num=0;
	int lyrs[7],sts[7],nbs[5];
	int rdm_nb;
	double rdm;

	double DB_temp1;
	double DB_temp2;

	rdm=RandRandN1()*diff_ttl;
	for(lyr=1;rdm>=0;lyr++)
		rdm-=p_diff[lyr];
	lyr--;
	rdm+=p_diff[lyr];
	for(st=0;rdm>=0;st++)
	{
		if(GET_SITE(lyr,st,0))
		{
			DB_temp1 = GET_NB_NUM(lyr,st);
			DB_temp2 = DIFF_RT[GET_NB_NUM(lyr,st)];
			rdm-=DIFF_RT[GET_NB_NUM(lyr,st)];
		}
	}
	st--;

	remove_particle(lyr,st);

	get_relative_sites(lyr, st, lyrs, sts);
	for(n=1;n<7;n++)
	{
		if((!GET_SITE(lyrs[n],sts[n],0))&&(GET_NB_NUM(lyrs[n],sts[n])>=1))
		{
			nbs[nb_num]=n;
			nb_num++;
		}
	}
	rdm_nb=int(RandRand()*nb_num);
	add_relax(lyrs[nbs[rdm_nb]],sts[nbs[rdm_nb]]);

	return 1;
}
//*******************************************************************


//*******************************************************************
// Utilities/Random generator
//*******************************************************************
double ranO(int *idum)
/*Returns a uniform random deviate between 0.0 and 1.0 using rand().*/
/*0.0 is available but 1.0 is not*/
/*Set idum to any negative value to initialize or reinitialize the sequence*/
{
	static double y, maxran, v[131];
	double dum;
	static int iff=0;
	int j;

	if (*idum < 0 || iff == 0)
	{
		iff = 1;
		maxran = RAND_MAX+0.0;
		srand((unsigned)time(NULL));
		*idum = 1;
		for (j=0; j<=130; j++)
			dum=rand()+0.0;
		for (j=0; j<=130; j++)
			v[j]=rand()+0.0;
		y=rand()+0.0;
	}
	j= (int)(130.0*y/maxran);
	if (j > 130 || j < 0)
		printf("RanO: This cannot happen.\n");
	y=v[j];
	v[j]=rand();
	return y/(maxran+1.0);
}
//*******************************************************************
double RandRand()
{
	double r_tot, r;
	int i, n;
	static int seed;

	n=3;

	r_tot=0.0;

	for (i=0; i<n; i++)
	{
		r=ranO(&seed);
		if (i > 0)
			r/=pow((RAND_MAX+1.0), (i+0.0));
		r_tot+=r;
	}
	
	return r_tot;
}
//*******************************************************************
double RandRandN1()
{
	double rdm;

	rdm=RandRand();
	while(rdm>=1){rdm=RandRand();}
	return rdm;
}
//*******************************************************************
double RandRandN0()
{
	double rdm;

	rdm=RandRand();
	while(rdm<=0){rdm=RandRand();}
	return rdm;
}
//*******************************************************************



//*******************************************************************
// Private functions/Simulation control
//*******************************************************************
double lattice_1DTRI::OneMCTime(double r_p)
{
	double rd;

	rd=RandRand();
	while ((rd <= 0) || (rd >= 1))
		rd=RandRand();

	return -log(rd)/r_p;
}
//*******************************************************************


//*******************************************************************
//
//	UNUSED FUNCTIONS
//
//*******************************************************************
/*
//*****************************************************************************
int lattice_1DTRI::test_lattice_cons()
{
	int i,j,n;
	int lyrs[7], sts[7];
	int cons=1; // all site values on the lattice are consistent

	for(j=1;j<=MLyr;j++)
	{
		for(i=0;i<LatticeSize;i++)
		{
			get_relative_sites(j, i, lyrs, sts);
			for(n=1;n<7;n++)
			{
				cons&=(GET_SITE(lyrs[0],sts[0],n)==GET_SITE(lyrs[n],sts[n],0));
			}
			if(GET_SITE(j,i,0)&&(GET_NB_NUM(j,i)==0))
				printf("Singlton appears!\n");
			if(GET_SITE(j,i,0)&&(GET_SITE(lyrs[1],sts[1],0)==0)&&(GET_SITE(lyrs[2],sts[2],0)==0)&&(GET_SITE(lyrs[3],sts[3],0)==0)&&(GET_SITE(lyrs[4],sts[4],0)==0)&&(GET_SITE(lyrs[5],sts[5],0)==0)&&(GET_SITE(lyrs[6],sts[6],0)==0))
				printf("Singlton appears!\n");
		}
	}

	if(cons)
	{
		printf("Consistent!\n");
	}
	else
	{
		printf("NOT Consistent!\n");
	}

	return cons;
}
*/
//*****************************************************************************

//*****************************************************************************
/*void lattice_1DTRI::save_lattice(int repeat, int Recorder_tIDX)
{
	char str[25],filename[80];  
	FILE *fp;
	int lyr,st;

	strcpy(filename, filename_l);
	strcat(filename, "_r");
	_itoa(repeat,str,10);
	strcat(filename, str);
	strcat(filename, "_t");
	_itoa(Recorder_tIDX,str,10);
	strcat(filename, str);
	initiate_file(filename);

	fp=fopen(filename, "a+");
	for(lyr=0;lyr<=TopLyr;lyr++)
	{
		for(st=0;st<LatticeSize;st++)
		{
			fprintf(fp, "%d ", GET_SITE(lyr,st,0));
		}
		fprintf(fp, "\n");
	}
	fclose(fp);

	return;
}*/

//*****************************************************************************
double lattice_1DTRI::calc_SOR()
{
	int lyr,st;
	int g_occu_sum=0,s_occu_sum=0;
	int top=0;
	int ht_seg;

	for(st=0;st<LatticeSize;st++)
	{
		for(hts[st]=MLyr;!GET_SITE(hts[st],st,0);hts[st]--) {}
		if(hts[st]>top)
		{
			top=hts[st];
		}
	}

	// Update TopLyr to top
	TopLyr=top;

	ht_seg=top>SURF?top:SURF;
	ht_seg-=SURF;

	for(st=0;st<LatticeSize;st++)
	{
		occu[st]=0;
		for(lyr=hts[st];lyr>ht_seg;lyr--)
		{
			occu[st]+=GET_SITE(lyr,st,0);
		}
		s_occu_sum+=occu[st];

		for(lyr=ht_seg;lyr>0;lyr--)
		{
			occu[st]+=GET_SITE(lyr,st,0);
		}
		g_occu_sum+=occu[st];
	}
	
	// to calculate the surface density, the surface layers are complemented to SURF layers
	// no need for complementation if the total layers are enough (top>=SURF)
	if(top<SURF)
	{
		s_occu_sum+=LatticeSize*(SURF-top);
	}
	LState.pSOR=(s_occu_sum+0.0)/(LatticeSize+0.0)/(SURF+0.0);
	
	if(top==0)
	{
		LState.SOR=0.0;
	}
	else
	{
		LState.SOR=(g_occu_sum+0.0)/(LatticeSize+0.0)/(top+0.0);
	}

	// Do NOT accumulate record.H and record.SOR here, since kMC_calc_stats() may be called by controller for actuation calculation

	return LState.SOR;
}
//*****************************************************************************
void lattice_1DTRI::calc_rough()
{
	int st;
	int init_lyr;
	double *db_hts; // used to save the surface height


	int i,j;

	db_hts = new double[2*LatticeSize];

	if(MLyr&0x1)
	{
		init_lyr=MLyr;
	}
	else
	{
		init_lyr=MLyr+1;
	}

	// find top position for each column and semi-column
	for(st=0;st<LatticeSize;st++)
	{
		for(hts[st]=init_lyr;!GET_SITE(hts[st],st,0);hts[st]-=2) 
		{
			if(hts[st]==1)
			{
				hts[st]=0;
				break;
			}
		}
		for(hts_half[st]=init_lyr-1;!GET_SITE(hts_half[st],st,0);hts_half[st]-=2) 
		{
			if(hts_half[st]==1)
			{
				hts_half[st]=0;
				break;
			}
		}
	}
	
	for(st=0;st<LatticeSize;st++)
	{
		db_hts[st*2]=hts_half[st];
		db_hts[st*2+1]=hts[st];
		if((hts_half[st]>hts[st])&&(hts_half[NEXT_ST[st]]>hts[st]))
		{
			db_hts[st*2+1]=(hts_half[st]+hts_half[NEXT_ST[st]])/2.0;
		}
		else if((hts[PREV_ST[st]]>hts_half[st])&&(hts[st]>hts_half[st]))
		{
			db_hts[st*2]=(hts[PREV_ST[st]]+hts[st])/2.0;
		}
	}
	
	// Calculate thickness
	double ht_sum=0.0;
	for(st=0;st<2*LatticeSize;st++)
	{
		ht_sum+=db_hts[st];
	}

	LState.H=ht_sum/(2*LatticeSize+0.0);
	
	// Calculate roughness
	double r2_sum=0.0;
	for(st=0;st<2*LatticeSize;st++)
	{
		r2_sum+=(db_hts[st]-LState.H)*(db_hts[st]-LState.H);
	}

	LState.R2=r2_sum/(2*LatticeSize+0.0);

	// Calculate root mean square slope
	double m2_sum = 0.0;
	for(st=0;st<2*LatticeSize-1;st++)
	{
		m2_sum+=(db_hts[st+1]-db_hts[st])*(db_hts[st+1]-db_hts[st])*3;
	}
//  m_sum+=m_sum+(db_hts[1]-db_hts[2*L_X])*(db_hts[1]-db_hts[2*L_X])/4+(db_hts[0]-db_hts[2*L_X-1])*(db_hts[0]-db_hts[2*L_X-1])/4;
    m2_sum=m2_sum+(db_hts[0]-db_hts[2*LatticeSize-1])*(db_hts[0]-db_hts[2*LatticeSize-1])*3	;
	LState.M2=m2_sum/(2*LatticeSize+0.0);

	// Calculate alpha and beta from the height profile using:
	// alpha_n=int_{-pi}^{pi}{phi_n(x)*h(x)*dx}=INT_COL*sum_0^L{phi_n(i)*h(i)}
	// beta_n=int_{-pi}^{pi}{psi_n(x)*h(x)*dx}=INT_COL*sum_0^L{psi_n(i)*h(i)}
	for (i=0; i<MODE; i++)
	{
		LState.alphas[i]=0.0;
		LState.betas[i]=0.0;
		for (j=0; j<2*LatticeSize; j++)
		{
			LState.alphas[i]+=db_hts[j]*phi_nx[i*2*LatticeSize+j];
			LState.betas[i]+=db_hts[j]*psi_nx[i*2*LatticeSize+j];
		}
		LState.alphas[i]*=INT_COL;
		LState.betas[i]*=INT_COL;
	}
	
	delete [] db_hts;
	return;
}
//*****************************************************************************
// keep updating records until:
// 1. all records are all updated (kMC simulation may finish) or
// 2. all possible records are updated right before iCurTime
void lattice_1DTRI::update_records (double iCurTime)
{
	int m;
	int i;
	dep_time=iCurTime;
	try
	{
		while(record.rec_idx<record.rec_num && record.pRecords[record.rec_idx].time<=iCurTime)
		{
			// calculate LState.R2 and LState.SOR
			// 2009.08.23 Xinyu: calc_rough() first, then calc_SOR() because for debug purpose, calc_SOR() will LState.H, which is calculated in calc_rough().
			calc_rough(); 
			calc_SOR();
			// Accumulate here
			record.pRecords[record.rec_idx].T			+=Temperature;
			record.pRecords[record.rec_idx].w			+=Dep_rate;
			record.pRecords[record.rec_idx].H_mean		+=LState.H;
			record.pRecords[record.rec_idx].H_var		+=LState.H*LState.H;
			record.pRecords[record.rec_idx].R2_mean	    +=LState.R2;
			record.pRecords[record.rec_idx].R2_var		+=LState.R2*LState.R2;
			record.pRecords[record.rec_idx].SOR_mean	+=LState.SOR;
			record.pRecords[record.rec_idx].SOR_var		+=LState.SOR*LState.SOR;
			record.pRecords[record.rec_idx].pSOR_mean	+=LState.pSOR;
			record.pRecords[record.rec_idx].pSOR_var	+=LState.pSOR*LState.pSOR;
			record.pRecords[record.rec_idx].M2_mean		+=LState.M2;
			record.pRecords[record.rec_idx].M2_var		+=LState.M2*LState.M2;

			if(Detail==1)
			{
				fd_time << record.pRecords[record.rec_idx].time << "\n";
				fd_T    << Temperature << "\n";
				fd_W    << Dep_rate << "\n";
				fd_SOR  << LState.SOR << "\n";
				fd_H    << LState.H << "\n";
				fd_R2   << LState.R2 << "\n";
				fd_M2   << LState.M2 << "\n";
				//for(i=0;i<MODE-1;i++)
				//{
					//fd_alpha << setw(15) << LState.alphas[i] << "\t";
					//fd_beta  << setw(15) << LState.betas[i] << "\t";
				//}
				//fd_alpha << setw(15) << LState.alphas[MODE-1] << "\n";
				//fd_beta  << setw(15) << LState.betas[MODE-1] << "\n";
			}
			record.rec_idx++;
		}
	}
	catch(const char *s)
	{
		cerr << "\n" << s;
		cerr << "\nError in update_records(): record.rec_idx = " << record.rec_idx;
	}

	return;
}
//*****************************************************************************
void lattice_1DTRI::avg_stats()
{
	int i;

	for (i=0;i<record.rec_num;i++)
	{
		record.pRecords[i].T			/=Repeated;
		record.pRecords[i].w			/=Repeated;		
		record.pRecords[i].H_mean		/=Repeated;
		record.pRecords[i].H_var		=record.pRecords[i].H_var/Repeated-record.pRecords[i].H_mean*record.pRecords[i].H_mean;
		record.pRecords[i].R2_mean		/=Repeated;
		record.pRecords[i].R2_var		=record.pRecords[i].R2_var/Repeated-record.pRecords[i].R2_mean*record.pRecords[i].R2_mean;
		record.pRecords[i].SOR_mean		/=Repeated;
		record.pRecords[i].SOR_var		=record.pRecords[i].SOR_var/Repeated-record.pRecords[i].SOR_mean*record.pRecords[i].SOR_mean;
		record.pRecords[i].pSOR_mean	/=Repeated;
		record.pRecords[i].pSOR_var		= record.pRecords[i].pSOR_var/Repeated-record.pRecords[i].pSOR_mean*record.pRecords[i].pSOR_mean;
		record.pRecords[i].M2_mean	/=Repeated;
		record.pRecords[i].M2_var	=record.pRecords[i].M2_var/Repeated-record.pRecords[i].M2_mean*record.pRecords[i].M2_mean;
	}
	return;
}
//*****************************************************************************
void lattice_1DTRI::getSnapshot(const string filename)
{
	ofstream of_Snapshot(filename.c_str());

	of_Snapshot << "LatticeType = 1D Triangular\n";
	of_Snapshot << "LatticeSize = " << LatticeSize << "\n";
	of_Snapshot << "TopLayer    = " << TopLyr << "\n";
	int i,j;
	// To save output 
	for(i = 0;i<TopLyr;i++)
	{
		for(j = 0;j<LatticeSize;j++)
		{
			if(!(GET_SITE(i,j,0))) {of_Snapshot << i << "\t" << j << "\n";}
		}
	}
	
	of_Snapshot.close();
}